# data vizlibrary(ggplot2)library(ggh4x)library(ggOceanMaps)library(patchwork)library(viridis)library(ggpubr)library(grid)library(viridis)library(ggnewscale)library(gridExtra)library(png)# tablelibrary(gt)library(gtable)# maplibrary(ggmap)library(ggsn)library(sf)library(sp)library(smoothr)# kernel densitylibrary(eks)library(ggdensity)# statlibrary(Hmisc)library(circular)library(CircStats)# charlibrary(stringr)# solar angle and timezonelibrary(suncalc)library(lutz)# data wranglinglibrary(tidyr)library(dplyr)library(data.table)library(purrr)library(forcats)library(lubridate)
2 Setting up custom function
2.1matlab_to_posix
Code
# function to convert matlab date into R posixmatlab_to_posix <-function(x, timez ="UTC") { days <- x -719529# 719529 = days from 1-1-0000 to 1-1-1970 secs <- days *86400# 86400 seconds in a dayreturn(as.POSIXct(strftime(as.POSIXct(secs,origin ="1970-1-1",tz ="UTC" ),format ="%Y-%m-%d %H:%M:%S",tz ="UTC",usetz =FALSE ),tz = timez ))}
2.2windrose
Code
windrose <-function(data_to_plot,grid =NULL,set_title =NULL,legend_position ="none",facet_wrap = F,facet_grid = F,hour_sunset =20,hour_sunrise =7) {# this code comes from Roxanne uniqhours <-1:24* (360/24)# trick to align hours om the graph data_to_plot <-rbind(data_to_plot[2:nrow(data_to_plot), ], data_to_plot[1, ])for (i in1:24) {# turn hours to radiansif (i ==1) { temp <-rep(uniqhours[i], data_to_plot$nb_ind_hour[i]) day_night <-rep("night", data_to_plot$nb_ind_hour[i]) } else { temp <-c(temp, rep(uniqhours[i], data_to_plot$nb_ind_hour[i])) day_night <-c(day_night, rep(if_else(between(i, hour_sunrise, hour_sunset), "day", "night"), data_to_plot$nb_ind_hour[i] )) } } data2 <-data.frame(direction = temp) deg <-15# choose bin size (degrees/bin) dir.breaks <-seq(0- (deg /2), 360+ (deg /2), deg) # define the range of each bin# assign each direction to a bin range dir.binned <-cut(data2$direction,breaks = dir.breaks,ordered_result =TRUE )# generate pretty lables dir.labels <-as.character(c(seq(0, 360- deg, by = deg), 0))# replace ranges with pretty bin labelslevels(dir.binned) <- dir.labels# Assign bin names to the original data set data2$dir.binned <- dir.binned# add origin if anyif (facet_wrap =="origin") { data2$origin <-unique(data_to_plot$origin) } elseif (facet_wrap =="trip") { data2$trip <-unique(data_to_plot$trip) }# add trip and origin if anyif (facet_grid) { data2$origin <-unique(data_to_plot$origin) data2$trip <-unique(data_to_plot$trip) }# set up max value maxvalue <-15# initialise the plot plt.dirrose_2 <-ggplot()# check if gridif (!is.null(grid)) { plt.dirrose_2 <- plt.dirrose_2 +geom_hline(yintercept = grid,colour ="grey20",linewidth = .2 ) } plt.dirrose_2 <- plt.dirrose_2 +geom_vline(xintercept =c(seq(1, 24, 2)),colour ="grey30",linewidth =0.2 ) +# 24 vertical lines at center of the 30? ranges.geom_hline(yintercept = maxvalue,colour ="black",linewidth = .5 ) +# Darker horizontal line as the top border (max).# On top of everything we place the histogram bars.geom_bar(data = data2,aes(x = dir.binned, fill = day_night),width =1,colour ="black",linewidth =0.3 ) +# geom_bar(data = data2, aes(x = dir.binned2), width = 1, colour="black", size = 0.3,fill="salmon",alpha=0.9) +scale_x_discrete(drop =FALSE,labels =c(0, "", 2, "", 4, "", 6, "", 8, "", 10, "", 12, "", 14, "", 16, "", 18, "", 20, "", 22, "" ) ) +scale_fill_manual(values =c("white", "darkgrey"), name ="Time of day") +labs(x ="Time (hours)", y ="Count", title = set_title) +# somehow was neededtheme(panel.grid.major =element_blank()) +coord_polar(start =-(deg /2) * (pi /180))# if facetif (facet_wrap =="origin") { plt.dirrose_2 <- plt.dirrose_2 +facet_wrap2(. ~ origin) } elseif (facet_wrap =="trip") { plt.dirrose_2 <- plt.dirrose_2 +facet_wrap2(trip ~ .,strip.position ="right" ) }# if facetif (facet_grid) { plt.dirrose_2 <- plt.dirrose_2 +facet_grid2(trip ~ origin) }# Wraps the histogram into a windrose plt.dirrose_2 <- plt.dirrose_2 +theme_bw() +theme(legend.position = legend_position,panel.grid.minor =element_blank(),panel.grid.major =element_blank(),panel.background =element_blank(),legend.key.size =unit(0.5, "cm") )# returnreturn(plt.dirrose_2) }
3 Import data
Let’s load data_dive, i.e. the output of data_wrangling.qmd, and filter only on animals leaving from Ano Nuevo.
Code
# import the data# data_dive <- readRDS("../export/data_dive_rev_3_inter.rds")data_dive <-readRDS("../export/data_dive.rds")# filter on seals departing from ANNUdata_dive <- data_dive %>%filter(DepartLoc =="ANNU")# remove id == `2015006` (no benthic dives)data_dive <- data_dive %>%filter(!id %in%c("2015006"))
We need to define a custom function to calculate the summarize version of each column differently, i.e. weighted mean and sum.
Code
compute_table_1_a <-function(x) {# if the sum of the column is below 0.1if (sum(x) < (0.1*2)) {# it's the proportion of dives that are benthic# so we weighted average these by taking into account all dives y <-weighted.mean(x, c( data_dive %>%filter(trip =="Post Breeding") %>%nrow(), data_dive %>%filter(trip =="Post Molting") %>%nrow() )) y <-paste0(round(y *100, 1), "%")return(y)# otherwise if below 1 } elseif (sum(x) < (1*2)) {# it's the proportion of benthic dives within predator areas# so we weighted average these by taking into account only dives with location data y <-weighted.mean(x, c( data_dive %>%filter(trip =="Post Breeding"&!is.na(Lat)) %>%nrow(), data_dive %>%filter(trip =="Post Molting"&!is.na(Lat)) %>%nrow() )) y <-paste0(round(y *100, 1), "%")return(y)# otherwise we just sum the column } else { y <-prettyNum(sum(x), big.mark =",")return(y) }}
# custom function to summarize each column differentlycompute_table_1_b <-function(x) {# if the column is numericif (all(is.numeric(x))) {# and the sum is below 100if (sum(x) <100) {# it is time# so we "circular" weighted average them y <- psych::circadian.mean(c(rep(x[1], data_dive %>%filter(trip =="Post Breeding"&!is.na(Lat)) %>%pull(id) %>%n_distinct()),rep(x[2], data_dive %>%filter(trip =="Post Molting"&!is.na(Lat)) %>%pull(id) %>%n_distinct()) ),na.rm = T ) y <-round(y, 1) } else {# otherwise it is just number# so we "just" weighted average them y <-weighted.mean(x, c( data_dive %>%filter(trip =="Post Breeding") %>%pull(id) %>%n_distinct(), data_dive %>%filter(trip =="Post Molting") %>%pull(id) %>%n_distinct() )) y <-round(y, 1) } } else { y <-"\u00b1" }return(y)}
Summary of dataset characteristics for Post Breeding and Post Molting
foraging trips in northern elephant seals.
# seals
All dives
Benthic dives
% Benthic dives per area
Total
w/ loc
w/o loc
Total
Proportion
Shark
Orca
Total
Overlap
Post Breeding
244
1,090,092
972,017
118,075
33,982
3.1%
67.5%
83.0%
83.2%
67.2%
Post Molting
233
2,597,260
1,980,478
616,782
51,748
2.0%
70.4%
87.2%
87.4%
70.2%
OVERALL
477
3,687,352
2,952,495
734,857
85,730
2.3%
69.4%
85.8%
86%
69.2%
Code
# display second tabletable_1_b
Table 2:
Individual-level summary of dataset characteristics for Post Breeding
and Post Molting foraging trips in northern elephant seals
# dives
# benthic dives
Hour departure
Hour arrival
Mean
±
SD
Mean
±
SD
Mean
±
SD
Mean
±
SD
Post Breeding
4467.6
±
964.6
139.3
±
230.7
20.6
±
1.7
9.3
±
2.8
Post Molting
11147.0
±
2093.1
222.1
±
384.3
19.3
±
1.7
13.7
±
1.9
OVERALL
7730.3
±
1515.8
179.7
±
305.7
20
±
1.7
11.2
±
2.4
Notes
Table 1: Summary of dataset characteristics for Post Breeding and Post Molting foraging trips in northern elephant seals.
Table 2: Individual-level summary of dataset characteristics for Post Breeding and Post Molting foraging trips in northern elephant seals
Code
# get the hours of departure and arrivalhour_data <- data_dive %>%group_by(trip, id) %>%# create a column to say when start and end each trip by idmutate(is_arrival =if_else(DiveNumber ==max(DiveNumber), T, F),is_departure =if_else(DiveNumber ==min(DiveNumber), T, F) ) %>%# testsummarise(hour_day_departure =hour(date_tz[is_departure == T]) +minute(date_tz[is_departure == T]) /60+second(date_tz[is_departure == T]) / (60*60),hour_day_arrival =hour(date_tz[is_arrival == T]) +minute(date_tz[is_arrival == T]) /60+second(date_tz[is_arrival == T]) / (60*60),.groups ="drop" ) %>%filter(!is.na(hour_day_departure) &!is.na(hour_day_arrival))# compare the circular mean hour of departure between two types of tripwatson.williams.test(list(pb = hour_data %>%filter(trip =="Post Breeding") %>%pull(hour_day_departure) %>%circular(., units ="hours"),pm = hour_data %>%filter(trip =="Post Molting") %>%pull(hour_day_departure) %>%circular(., units ="hours") ))
Watson-Williams test for homogeneity of means
data: pb and pm
F = 2.8203, df1 = 1, df2 = 214, p-value = 0.09454
sample estimates:
Circular Data:
Type = angles
Units = hours
Template = none
Modulo = asis
Zero = 0
Rotation = counter
mean of pb mean of pm
-2.819599 -4.515445
Code
watson.wheeler.test(list(pb = hour_data %>%filter(trip =="Post Breeding") %>%pull(hour_day_departure) %>%circular(., units ="hours"),pm = hour_data %>%filter(trip =="Post Molting") %>%pull(hour_day_departure) %>%circular(., units ="hours") ))
Watson-Wheeler test for homogeneity of angles
data: pb and pm
W = 1.1476, df = 2, p-value = 0.5634
Individual are not chosen completely randomly, we picked randomly individuals that performed ~180 dives per 36h.
Code
# to choose an individual ~180 dives for 36 hoursdata_dive %>%# then by individualgroup_by(id) %>%# calculatemutate( # the difference time with the first divediff_start =abs(as.numeric(difftime(first(DateTime), DateTime,units ="days" ))),# the difference time with the last divediff_end =abs(as.numeric(difftime(last(DateTime), DateTime,units ="days" ))) ) %>%filter(diff_start <= nb_days | diff_end <= nb_days) %>%mutate(phase =factor(if_else(diff_start <= nb_days, "Departure", "Arrival"),levels =c("Departure", "Arrival") )) %>%group_by(id, trip) %>%summarise(nb_dives =length(DiveNumber)) %>%arrange(nb_dives) %>%View()
Code
# import the csv exporttdr_data <-lapply(list.files(path ="../data/",pattern ="*_iknos_raw_data.csv",full.names = T ),function(x) {# extract id id_name <-last(strsplit(strsplit(x, split ="_")[[1]][1],split ="/" )[[1]])# load file dt <- readr::read_csv(x)# add id dt %>%mutate(id =as.numeric(id_name)) }) %>%bind_rows()# choose the lagnb_days <-1.5# compute the datadata_plot <- tdr_data %>%# convert time into date, a Posixmutate(DateTime =matlab_to_posix(time),date =as.Date(DateTime) ) %>%# add beginning and end date of the tripmerge( ., data_dive %>%group_by(id, trip) %>%summarise(date_low =min(DateTime),date_high =max(DateTime),.groups ="drop" ),by ="id",all.x =TRUE ) %>%# identify what to keepmutate(to_keep =if_else(between(DateTime, date_low, date_high), T, F)) %>%# keep only between the beginning and the endfilter(to_keep) %>%# sort to make sure the difftime are being done correctlyarrange(id, time) %>%# merge to get lat and lonmerge( ., data_dive %>%group_by(id, date) %>%summarise(lat =median(Lat, na.rm = T),lon =median(Lon, na.rm = T),.groups ="drop" ),by =c("id", "date"),all.x =TRUE ) %>%# then by individualgroup_by(id) %>%# calculatemutate( # the difference time with the first divediff_start =abs(as.numeric(difftime(first(DateTime), DateTime,units ="days" ))),# the difference time with the last divediff_end =abs(as.numeric(difftime(last(DateTime), DateTime,units ="days" ))) ) %>%filter(diff_start <= nb_days | diff_end <= nb_days) %>%mutate(phase =factor(if_else(diff_start <= nb_days, "Departure", "Arrival"),levels =c("Departure", "Arrival") ))# add bathysetDT(data_plot)setDT(data_dive)setkeyv(data_plot, c("id", "DateTime"))setkeyv(data_dive, c("id", "DateTime"))data_plot <-copy(data_dive) %>% .[, .(id, DateTime, bathy, DiveTypeName)] %>% .[data_plot, roll = T]# get sunrise sunset informationsun_data <- data_plot %>%select(date, lat, lon, id, phase) %>%unique() %>%ungroup() %>%mutate(getSunlightTimes(data = ., keep =c("sunrise", "sunset"))) %>%# fuck it, getSunlightTimes gives the sunset of next day... so we subtract 1mutate(sunset =if_else( # if sunset is the next dayabs(as.numeric(difftime(as.POSIXct(date), sunset, units ="days") )) >1,# subtract 1 sunset - (1*60*60*24),# otherwise leave it that way sunset ))# trim if night beging before the trip or end after they arrivesun_data <-merge( sun_data %>%select(date, lat, lon, id, phase, sunrise, sunset), data_plot %>%group_by(id, phase, trip) %>%summarise(date_low =min(DateTime),date_high =max(DateTime),.groups ="drop" ),by =c("id", "phase"),all.x =TRUE, ) %>%mutate(# make sure sunset is not outside the 36h periodcorrect_sunset =if_else(# if notbetween(sunset, date_low, date_high), sunset,# if so, replace by one extremumif_else(sunset < date_low, date_low, date_high) ),# make sure sunrise is not outside the 36h periodcorrect_sunrise =if_else(# if notbetween(sunrise, date_low, date_high), sunrise,# if so, replace by one extremumif_else(sunrise < date_low, date_low, date_high) ) ) %>%# shift everything to start 2000-01-01 00:00:00mutate(correct_sunset_shift = correct_sunset - ( date_low -as.POSIXct("2000-01-01 00:00:00", tz ="UTC") ),correct_sunrise_shift = correct_sunrise - ( date_low -as.POSIXct("2000-01-01 00:00:00", tz ="UTC") ) )# plotdeparture_dp <- data_plot %>%filter(phase =="Departure") %>%group_by(id, phase) %>%# make the date starts on 2000-01-01 00:00:00mutate(DateTime = DateTime - (first(DateTime) -as.POSIXct("2000-01-01 00:00:00", tz ="UTC") )) %>%ungroup() %>%mutate(`Dive Type`=if_else(DiveTypeName =="Benthic", "Benthic", "Other")) %>%ggplot(aes(x = DateTime, y = depth, col =`Dive Type`, group = id)) +scale_color_manual(values =c("Benthic"="#008c49", "Other"="grey10")) +geom_rect(data = sun_data %>%filter(phase =="Departure"),aes(xmin = correct_sunset_shift,xmax = correct_sunrise_shift,ymin =-Inf,ymax =+Inf ),alpha =0.4,inherit.aes =FALSE ) +# geom_path(aes(x = DateTime, y = abs(bathy))) +geom_path() +scale_x_datetime(# date_breaks = "6 hours",# date_labels = "+ %Hh",breaks =seq.POSIXt(as.POSIXct("2000-01-01 06:00:00", tz ="UTC"),as.POSIXct("2000-01-02 06:00:00", tz ="UTC"),by ="6 hour" ),labels =c("+6h", "+12h", "+18h", "+24h", "+30h"),position ="top" ) +scale_y_reverse(limits =rev(range(data_plot$depth))) +labs(x ="Time", y ="Depth (m)") +# facet_nested(trip + id ~ phase,facet_nested(id ~ phase,scales ="free_x",# independent = "x" ) +theme(strip.placement ="outside",strip.text.y =element_blank(),legend.position ="top",panel.grid.minor =element_blank(),panel.border =element_blank(),panel.background =element_blank(),panel.grid.major =element_line(color ="grey90"),axis.line.x =element_line(color ="black",# arrow = arrow(length = unit(0.2, "lines"), type = "closed") ),axis.line.y =element_line(color ="black",arrow =arrow(length =unit(0.2, "lines"),type ="closed",ends ="first" ) ),axis.title.x =element_blank() )arrival_dp <- data_plot %>%filter(phase =="Arrival") %>%group_by(id, phase) %>%# make the date starts on 2000-01-01 00:00:00mutate(DateTime = DateTime - (first(DateTime) -as.POSIXct("2000-01-01 00:00:00", tz ="UTC") )) %>%ungroup() %>%mutate(`Dive Type`=if_else(DiveTypeName =="Benthic", "Benthic", "Other")) %>%ggplot(aes(x = DateTime, y = depth, col =`Dive Type`, group = id)) +scale_color_manual(values =c("Benthic"="#008c49", "Other"="grey10")) +geom_rect(data = sun_data %>%filter(phase =="Arrival"),aes(xmin = correct_sunset_shift,xmax = correct_sunrise_shift,ymin =-Inf,ymax =+Inf ),alpha =0.4,inherit.aes =FALSE ) +# geom_path(aes(x = DateTime, y = abs(bathy))) +geom_path() +scale_x_datetime(# date_breaks = "6 hours",# date_labels = "+ %Hh",breaks =seq.POSIXt(as.POSIXct("2000-01-01 06:00:00", tz ="UTC"),as.POSIXct("2000-01-02 06:00:00", tz ="UTC"),by ="6 hour" ),labels =c("-30h", "-24h", "-18h", "-12h", "-6h"),position ="top" ) +scale_y_reverse(limits =rev(range(data_plot$depth))) +labs(x ="Time", y ="Depth (m)") +# facet_nested(trip + id ~ phase,facet_nested(id ~ phase,scales ="free_x",# independent = "x" ) +theme(strip.placement ="outside",legend.position ="top",panel.grid.minor =element_blank(),panel.border =element_blank(),panel.background =element_blank(),panel.grid.major =element_line(color ="grey90"),axis.line.x =element_line(color ="black",arrow =arrow(length =unit(0.2, "lines"), type ="closed") ),axis.line.y =element_blank(),axis.ticks.y =element_blank(),axis.text.y =element_blank(),axis.title.y =element_blank(),axis.title.x =element_blank() )fig_label_x <-ggplot(data.frame(l ="Time", x =1, y =1)) +geom_text(aes(x, y, label = l)) +theme_void() +coord_cartesian(clip ="off")
Figure 1: Snapshot of departure and arrival dive profiles for six seals. Departure and arrival were constrained to 36 hours. Shading indicates night and coloration indicates dive type.
Notes
TODO: Add three dots between Departure and Arrival
4.3 Figure 2
Let’s create a table specific for this figure that only contains for each individual:
the first dive
the last dive
all benthic dives
Code
# let's add a column with the local_timedata_windrose <- data_dive %>%# first we need to identify the start and the end of each trip by idgroup_by(id, trip) %>%mutate(is_arrival =if_else(DiveNumber ==max(DiveNumber), T, F),is_departure =if_else(DiveNumber ==min(DiveNumber), T, F) ) %>%# keep only the data of the first, last divesfilter(is_arrival == T | is_departure == T) %>%# ungroupungroup() %>%# then, since we're looking at local time, we keep data with locationfilter(!is.na(Lat)) %>%# change colnames to make getSunlightTimes worksmutate(date =as.Date(date_tz),lat = Lat,lon = Lon ) %>%# calculate hour of sunset and hour of sunrisemutate(getSunlightTimes(data = .,keep =c("sunrise", "sunset"),tz =tz_lookup_coords(lat = Lat,lon = Lon,method ="accurate" ) )) %>%# extract hours of sunset an sunrisemutate(h_sunrise =round(hour(sunrise) +minute(sunrise) /60+second(sunrise) / (60*60)),h_sunset =round(hour(sunset) +minute(sunset) /60+second(sunset) / (60*60)) )
Figure 2: Circular histogram plots display the local timing (in hours) for the first and last recorded dives (respectively departure and arrival) of adult female northern elephant seals (n = 403) during their Post Breeding and Post Molting foraging trips.
Note
n = 361, this is the number of individuals having location data (allowing to calculate the local time) for their first and/or their last dives
# functio to convert time in radian (probably already exist in some package...)hour_time_2_radian <-function(x, high =24, low =0) {# https://pingouin-stats.org/build/html/generated/pingouin.convert_angles.html#pingouin.convert_angles ptp <- high - low rad <- x * (2* pi) / ptp rad <- (rad + pi) %% (2* pi) - pireturn(rad)}
(Not meant to be in the paper) Average hour of departure and arrival
for post breeding and post molting trip with the result of the
associated Rayleigh test.
Figure 3: Spatial overlap of benthic dives performed by northern elephant seals (n = 401), and the shark (yellow polygon) and orca (red polygon) distribution area along the west coast (adapted from Jorgenson, et al. 2019) with a kernel density estimation representing the occurrence of all the benthic dives.
Figure 4: Same as the previous one, but with the benthic dives (black dots). We can see that some benthic happen outside orca areas…
Notes
?
n = 392, this is the number of individuals having at least one benthic with one location associated with (filter(!is.na(Lat) & DiveTypeName == "Benthic"))
4.5 Figure 4
Code
# https://mycolor.space/?hex=%23008C49&sub=1# build palette colorcols_ind_1 <-c("#008c49", "#00c9d3", "#bffcfd", "#458081")# test for color blind# => https://davidmathlogic.com/colorblind/#%23008C49-%2300C9D3-%23BFFCFD-%23458081cols_ind_2 <-c("#008c49", "#002e00", "#00c9d3", "#deb887")# test for color blind# => https://davidmathlogic.com/colorblind/#%23008C49-%23002E00-%2300C9D3-%23DEB887
Code
# main plotbathy_prop <- data_dive %>%# keep only rows with location datafilter(!is.na(Lat)) %>%# keep only negative bathyfilter(bathy <0) %>%# and remove outliersfilter(bathy >-6000) %>%# positive bathymutate(bathy =abs(bathy)) %>%# create class of bathymetrymutate(bathy_class =cut( bathy,seq(0, 6000, 400),ordered_result = T,dig.lab =4 )) %>%# calculate by bath_class and animalgroup_by(bathy_class, DiveTypeName) %>%# the number of divessummarise(N =n(), .groups ="drop") %>%# and divide by the total number of dives per bathy_classgroup_by(bathy_class) %>%# to get the Proportion of different dive types per bathy_classmutate(Proportion = N /sum(N)) %>%# ungroup => not required but that let the dataset cleanungroup() %>%# then plotggplot(aes(x = bathy_class, y = Proportion)) +# the areageom_area(aes(fill = DiveTypeName, group = DiveTypeName)) +# orientation of x labelsguides(x =guide_axis(angle =45)) +# format y axisscale_y_continuous(labels =function(x) {paste0(x *100, "%") } ) +labs(x ="Bathymetry (m)",y ="Dive type proportion (%)" ) +# scale_fill_viridis(# "Dive Type",# option = "plasma",# discrete = T,# direction = -1# ) +scale_fill_manual("Dive Type",values = cols_ind_1 ) +# position legend at the toptheme(legend.position ="top",axis.title.x =element_text(margin =margin(t =15)),panel.grid.major =element_blank(),panel.grid.minor =element_blank(),panel.background =element_blank(),axis.line =element_line(arrow =arrow(angle =20,length =unit(.10, "inches"),type ="closed" )),plot.tag.position =c(-0.03, 1) )# # second plot# bathy_hist <- data_dive %>%# # keep only negative bathy# filter(bathy < 0) %>%# # and remove outliers# filter(bathy > -6000) %>%# # create class of bathymetry# mutate(bathy_class = fct_rev(cut(# bathy,# seq(-6000, 0, 400),# ordered_result = T,# dig.lab = 4# ))) %>%# # calculate by bath_class and animal# group_by(bathy_class, DiveTypeName) %>%# # the number of dives# summarise(N = n(), .groups = "drop") %>%# # ggplot# ggplot(aes(# x = bathy_class,# y = N,# fill = DiveTypeName,# width = 0.5# )) +# geom_bar(# stat = "identity",# position = "dodge"# ) +# scale_fill_viridis(# option = "plasma",# discrete = T,# direction = -1# ) +# theme_void() +# theme(legend.position = "top")## # the actual plot# bath_plot <- bathy_hist / bathy_prop + plot_layout(heights = c(1, 10))
Code
# main plotdist_coast_prop <- data_dive %>%# filter out animal without dist_coastfilter(dist_coast >0) %>%# remove outliersfilter(dist_coast *1000*111<=1900) %>%# create class of bathymetrymutate(dist_class =cut(# convert decimal degree/1000 dist_coast *1000*111,seq(0, 1900, 100),ordered_result = T,dig.lab =4 )) %>%# calculate by bath_class and animalgroup_by(dist_class, DiveTypeName) %>%# the number of divessummarise(N =n(), .groups ="drop") %>%# and divide by the total number of dives per dist_classgroup_by(dist_class) %>%# to get the Proportion of different dive types per dist_classmutate(Proportion = N /sum(N)) %>%# ungroup => not required but that let the dataset cleanungroup() %>%# then plotggplot(aes(x = dist_class, y = Proportion)) +# the areageom_area(aes(fill = DiveTypeName, group = DiveTypeName)) +# orientation of x labelsguides(x =guide_axis(angle =45)) +# format y axisscale_y_continuous(labels =function(x) {paste0(x *100, "%") } ) +labs(x ="Distance from the coast (km)",y ="Dive type proportion" ) +# scale_fill_viridis(# "Dive Type",# option = "plasma",# discrete = T,# direction = -1# ) +scale_fill_manual("Dive Type",values = cols_ind_1 ) +# position legend at the toptheme(legend.position ="none",axis.title.x =element_text(margin =margin(t =15)),panel.grid.major =element_blank(),panel.grid.minor =element_blank(),panel.background =element_blank(),axis.line =element_line(arrow =arrow(angle =20,length =unit(.10, "inches"),type ="closed" )),plot.tag.position =c(-0.03, 1) )# # the second plot# dist_coast_hist <- data_dive %>%# # filter out animal without dist_coast# filter(dist_coast > 0) %>%# # remove outliers# filter(dist_coast * 1000 * 111 <= 1900) %>%# # create class of bathymetry# mutate(dist_class = cut(# # convert decimal degree/1000# dist_coast * 1000 * 111,# seq(0, 2000, 100),# ordered_result = T,# dig.lab = 4# )) %>%# # calculate by bath_class and animal# group_by(dist_class, DiveTypeName) %>%# # the number of dives# summarise(N = n(), .groups = "drop") %>%# # ggplot# ggplot(aes(# x = dist_class,# y = N,# fill = DiveTypeName,# width = 0.5# )) +# geom_bar(# stat = "identity",# position = "dodge"# ) +# scale_fill_viridis(# option = "plasma",# discrete = T,# direction = -1# ) +# theme_void() +# theme(legend.position = "top")## the actual plot# dist_coast_plot <- dist_coast_hist / dist_coast_prop + plot_layout(heights = c(1, 10))
Code
# the main plotdist_dep_prop <- data_dive %>%# filter out animal without dist_coastfilter(dist_dep >0) %>%# remove outliersfilter(dist_dep /1000<4200) %>%# create class of bathymetrymutate(dist_class =cut( dist_dep /1000,seq(0, 4200, 300),ordered_result = T,dig.lab =4 )) %>%# calculate by bath_class and animalgroup_by(dist_class, DiveTypeName) %>%# the number of divessummarise(N =n(), .groups ="drop") %>%# and divide by the total number of dives per dist_classgroup_by(dist_class) %>%# to get the Proportion of different dive types per dist_classmutate(Proportion = N /sum(N)) %>%# ungroup => not required but that let the dataset cleanungroup() %>%# then plotggplot(aes(x = dist_class, y = Proportion)) +# the areageom_area(aes(fill = DiveTypeName, group = DiveTypeName)) +# orientation of x labelsguides(x =guide_axis(angle =45)) +# format y axisscale_y_continuous(labels =function(x) {paste0(x *100, "%") } ) +labs(x ="Distance from A<U+00F1>o Nuevo (km)",y ="Dive type proportion" ) +# scale_fill_viridis(# "Dive Type",# option = "plasma",# discrete = T,# direction = -1# ) +scale_fill_manual("Dive Type",values = cols_ind_1 ) +# position legend at the toptheme(legend.position ="none",axis.title.x =element_text(margin =margin(t =15)),panel.grid.major =element_blank(),panel.grid.minor =element_blank(),panel.background =element_blank(),axis.line =element_line(arrow =arrow(angle =20,length =unit(.10, "inches"),type ="closed" )),plot.tag.position =c(-0.03, 1) )# # the second plot# dist_dep_hist <- data_dive %>%# # filter out animal without dist_coast# filter(dist_coast > 0) %>%# # remove outliers# filter(dist_dep / 1000 < 4200) %>%# # create class of bathymetry# mutate(dist_class = cut(# dist_dep / 1000,# seq(0, 4200, 300),# ordered_result = T,# dig.lab = 4# )) %>%# # calculate by bath_class and animal# group_by(dist_class, DiveTypeName) %>%# # the number of dives# summarise(N = n(), .groups = "drop") %>%# # ggplot# ggplot(aes(# x = dist_class,# y = N,# fill = DiveTypeName,# width = 0.5# )) +# geom_bar(# stat = "identity",# position = "dodge"# ) +# scale_fill_viridis(# option = "plasma",# discrete = T,# direction = -1# ) +# theme_void() +# theme(legend.position = "top")## # the actual plot# dist_dep_plot <- dist_dep_hist / dist_dep_prop + plot_layout(heights = c(1, 10))
Figure 5: Proportion of each dive type performed by northern elephant seals in relation to (A) bathymetry, (B) distance from the coast, and (C) distance from departure (defined as the first recorded location).
Notes
n = 406, because we removed individuals without location data (filter(!is.na(Lat)))
4.6 Figure 5
Code
# the last plotprop_at_sea_dt <- data_dive %>%group_by(id) %>%# time difference between the first date and the others (telling R to take diff b/w first date and all other dates)mutate(nb_days_departure =trunc(as.numeric(difftime( date, first(date),units ="days" )))) %>%# group by day/dive_type/sealgroup_by(nb_days_departure, DiveTypeName, id) %>%# count number of divessummarise(nb_daily_dives_type =n(), .groups ="drop") %>%# group by day/sealgroup_by(nb_days_departure, id) %>%# count the total number of dives per day/sealmutate(nb_daily_dives =sum(nb_daily_dives_type)) %>%# order by seals/datearrange(id, nb_days_departure) %>%# perform our calculation per sealgroup_by(id) %>%# calculate of the Proportion of time since departuremutate(perc_time_at_sea =round( nb_days_departure /max(nb_days_departure) *100, 1 )) %>%# filter on benthic divesfilter(DiveTypeName =="Benthic") %>%# add class of Proportion of time at seamutate(day_class =cut(perc_time_at_sea, seq(0, 100, 5),include.lowest = T )) %>%# by class of Proportion of time at seagroup_by(day_class) %>%# calculate the proportion of daily benthic divessummarise(nb_benthic_class =sum(nb_daily_dives_type),nb_total_class =sum(nb_daily_dives),prop_class = nb_benthic_class / nb_total_class )# plotprop_at_sea <- prop_at_sea_dt %>%# intiate plotggplot( .,aes(x = day_class, y = prop_class) ) +geom_bar(stat ="identity", fill ="#008c49") +guides(x =guide_axis(angle =45)) +scale_y_continuous(limits =c(0, 1),# breaks=seq(0,0.6,0.1),labels =function(x) { scales::percent(abs(x), 1) } ) +annotate("text", x =5.3, y =0.4, label ="Departure") +annotate("text", x =15, y =0.43, label ="Arrival") +# geom_curve(# data = tibble(# x1 = c(4, 17),# x2 = c(1.5, 19.5),# y1 = c(0.4, 0.5),# y2 = c(0.31, 0.55)# ),# aes(# x = x1,# y = y1,# xend = x2,# yend = y2# ),# arrow = arrow(length = unit(0.08, "inch")),# linewidth = 0.5,# color = "gray20",# curvature = 0.3# ) +geom_curve(data =tibble(x1 =c(4),x2 =c(1),y1 =c(0.4),y2 =c(0.33) ),aes(x = x1,y = y1,xend = x2,yend = y2 ),arrow =arrow(length =unit(0.08, "inch")),linewidth =0.5,color ="gray20",curvature =0.3 ) +geom_curve(data =tibble(x1 =c(16),x2 =c(19.7),y1 =c(0.43),y2 =c(0.5) ),aes(x = x1,y = y1,xend = x2,yend = y2 ),arrow =arrow(length =unit(0.08, "inch")),linewidth =0.5,color ="gray20",curvature =-0.56 ) +theme(panel.grid.major =element_blank(),panel.grid.minor =element_blank(),panel.background =element_blank(),axis.line =element_line(arrow =arrow(angle =20,length =unit(.10, "inches"),type ="closed" )) ) +labs(x ="Time spent at sea (%)",y ="Proportion of benthic dives" )
Code
# displayprop_at_sea
Figure 6: Proportion of benthic dives performed during trip at sea expressed as a percentage of time spent at sea. Time spent at sea is classified into intervals of five percent for graphical representation.
Notes
n = 477 all individuals (since they all have at leat one benthic dives)
Test to compare if the proportion of benthic dives at arrival is higher than the one at departure:
\(H_{0}: p_{departure} \ge p_{arrival}\)
\(H_{1}: p_{departure} \lt p_{arrival}\)
Code
# test to compare the proportion between arrival and departurez_test <-prop.test(x =c( prop_at_sea_dt %>%filter(day_class ==first(day_class)) %>%pull(nb_benthic_class), prop_at_sea_dt %>%filter(day_class ==last(day_class)) %>%pull(nb_benthic_class) ),n =c( prop_at_sea_dt %>%filter(day_class ==first(day_class)) %>%pull(nb_total_class), prop_at_sea_dt %>%filter(day_class ==last(day_class)) %>%pull(nb_total_class) ),alternative ="less")z_test
2-sample test for equality of proportions with continuity correction
data: c(prop_at_sea_dt %>% filter(day_class == first(day_class)) %>% pull(nb_benthic_class), prop_at_sea_dt %>% filter(day_class == last(day_class)) %>% pull(nb_benthic_class)) out of c(prop_at_sea_dt %>% filter(day_class == first(day_class)) %>% pull(nb_total_class), prop_at_sea_dt %>% filter(day_class == last(day_class)) %>% pull(nb_total_class))
X-squared = 4305.3, df = 1, p-value < 2.2e-16
alternative hypothesis: less
95 percent confidence interval:
-1.0000000 -0.1808821
sample estimates:
prop 1 prop 2
0.2880853 0.4734710
Code
# extract the z statisticsqrt(z_test$statistic)
X-squared
65.61442
\(H_{0}\) is rejected, so we can say that the porportion of benthic dives at arrival is significantly superior than departure!
# data use to compute kernel density estimationdf_kernel_dens <- data_dive %>%# only with location datafilter(!is.na(Lat)) %>%# select only the required columnsselect(lat = Lat, long = Lon, id, DiveTypeName) %>%# create col to nicely display dive typemutate(origin =paste(DiveTypeName, "dives"))
Code
# transform data into sf objectdf_kernel_dens_sf <-st_as_sf( df_kernel_dens,coords =c("long", "lat"),crs =st_crs(4326))# make it's group_by origindf_kernel_dens_sf <-group_by(df_kernel_dens_sf, origin)# kernel density estimationdf_kernel_dens_sf_kde <-st_kde( df_kernel_dens_sf,# H = diag(c(# MASS::bandwidth.nrd(# sf::st_coordinates(df_kernel_dens_sf)[, 1]# ),# MASS::bandwidth.nrd(# sf::st_coordinates(df_kernel_dens_sf)[, 2]# )# ) / 4)^2)# https://github.com/r-spatial/sf/issues/1762sf::sf_use_s2(FALSE)# plotto_plot <- trip_zoom_out +# rename axislabs(x ="Longitude", y ="Latitude") +# new scalenew_scale_fill() +# kernelgeom_sf(data =st_get_contour(# geospatial kernel df_kernel_dens_sf_kde,# probabilitiescont =c(50, 80, 95, 99) ),# displayaes(fill =label_percent(contlabel)),alpha =0.7 ) +# same colour barscale_fill_viridis_d(option ="plasma") +# legendlabs(fill ="Probability contours") +# no display alphaguides(alpha ="none") +# facet by originfacet_wrap(. ~ origin)# reorder layersto_plot$layers <-if (with_bathy) { to_plot$layers[c(1, 2, 4, 3)]} else { to_plot$layers[c(1, 3, 2)]}# displayto_plot
Figure 7: Map of kernel density estimation representing the spatial distribution of Benthic, Drift, Forage, and Transit dives performed by northern elephant seals throughout their trip to sea. Probability contours represent different levels of probability density.
Notes
lmk!
4.7.3 Figure 2
Code
# the last plotprop_at_sea_trip <- data_dive %>%group_by(id, trip) %>%# time difference between the first date and the others (telling R to take diff b/w first date and all other dates)mutate(nb_days_departure =trunc(as.numeric(difftime( date, first(date),units ="days" )))) %>%# group by day/dive_type/sealgroup_by(nb_days_departure, DiveTypeName, id, trip) %>%# count number of divessummarise(nb_daily_dives_type =n(), .groups ="drop") %>%# group by day/sealgroup_by(nb_days_departure, id) %>%# count the total number of dives per day/sealmutate(nb_daily_dives =sum(nb_daily_dives_type)) %>%# order by seals/datearrange(id, nb_days_departure, trip) %>%# perform our calculation per sealgroup_by(id, trip) %>%# calculate of the Proportion of time since departuremutate(perc_time_at_sea =round( nb_days_departure /max(nb_days_departure) *100, 1 )) %>%# filter on benthic divesfilter(DiveTypeName =="Benthic") %>%# add class of Proportion of time at seamutate(day_class =cut(perc_time_at_sea, seq(0, 100, 5),include.lowest = T )) %>%# by class of Proportion of time at seagroup_by(day_class, trip) %>%# calculate the proportion of daily benthic divessummarise(nb_benthic_class =sum(nb_daily_dives_type),nb_total_class =sum(nb_daily_dives),prop_class = nb_benthic_class / nb_total_class ) %>%# plotggplot( .,aes(x = day_class, y = prop_class) ) +geom_bar(stat ="identity", fill ="#008c49") +guides(x =guide_axis(angle =45)) +scale_y_continuous(limits =c(0, 1),# breaks=seq(0,0.6,0.1),labels =function(x) { scales::percent(abs(x), 1) } ) +facet_wrap(. ~ trip) +# annotate("text", x = 5.3, y = 0.4, label = "Departure") +# annotate("text", x = 15, y = 0.5, label = "Arrival") +# geom_curve(# data = tibble(# x1 = c(4),# x2 = c(1),# y1 = c(0.4),# y2 = c(0.33)# ),# aes(# x = x1,# y = y1,# xend = x2,# yend = y2# ),# arrow = arrow(length = unit(0.08, "inch")),# linewidth = 0.5,# color = "gray20",# curvature = 0.3# ) +# geom_curve(# data = tibble(# x1 = c(16),# x2 = c(19.7),# y1 = c(0.5),# y2 = c(0.57)# ),# aes(# x = x1,# y = y1,# xend = x2,# yend = y2# ),# arrow = arrow(length = unit(0.08, "inch")),# linewidth = 0.5,# color = "gray20",# curvature = -0.56# ) +theme(panel.grid.major =element_blank(),panel.grid.minor =element_blank(),panel.background =element_blank(),axis.line =element_line(arrow =arrow(angle =20,length =unit(.10, "inches"),type ="closed" )) ) +labs(x ="Time spent at sea (%)",y ="Proportion of benthic dives" )
Code
# displayprop_at_sea_trip
Figure 8: Proportion of benthic dives performed during trip at sea expressed as a percentage of time spent at sea during post breeding and post molting trip. Time spent at sea is classified into intervals of five percent for graphical representation.